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From the striped coats of zebras to the ripples in windblown sand, the natural 

ov 

world abounds with locally banded patterns. Such patterns have been of great interest 
throughout history, and, in the last twenty years, scientists in a wide variety of fields 
Mh. have been studying the patterns formed in well-controlled experiments that yield 

enormous quantities of high-precision data. These experiments involving phenomena 
as diverse as chemical reactions in shallow layers,El thermal convection in horizontal 
fluid layersjj periodically shaken layers of sandJ3 and the growth of slime mold colonies! 
often display patterns that appear qualitatively similar.u Methods are needed to 
characterize in a reasonable amount of time the differences and similarities in patterns 
that develop in different systems, as well as in patterns formed in one system for 
different experimental conditions. In this Letter, we introduce a novel, fast method for 
determining local pattern properties such as wavenumber, orientation, and curvature 
as a function of position for locally striped patterns. 

Fig. m shows a locally striped pattern formed in a simulation of Rayleigh-Benard convection, 
in which a horizontal layer of fluid of lateral extent L is confined between two parallel plates 
with separation d and a temperature difference AT = Tbottom — ^top > 0. For AT larger than 
a critical value AT C , a convective instability develops. In most situations studied experimentally, 



this instability leads to the development of a striped convective roll pattern with a horizontal local 
wavenumber k ~ ir/d. For systems of large aspect ratio T = L/d (i.e., many convective rolls), the 
patterns formed are often spatially and temporally disordered but are still locally striped as seen in 
Fig. [l|. The development of disordered locally periodic patterns through instabilities is common to 
a wide variety of large aspect ratio experiments in biological, chemical, optical, and fluid systems. 

The use of local pattern properties to describe patterns similar to that in Fig. [I| has already 
been explored by several researchers. Heutmaker and GollubQ determined the wavevector field and 
statistics such as the "roll bending" and "roll obliqueness" in both stationary and time-dependent 
patterns in a circular Rayleigh-Benard convection cell. They studied whether these quantities could 
be used to understand the stability of convective patterns. More recently, Hu, et alQ computed 
local wavenumbers and curvatures in experimental pictures very similar to the one shown in Fig. [l|. 
Based on these local measurements, they proposed order parameters to describe transitions in 
spatiotemporal chaos in Rayleigh-Benard convection. In a numerical study of a model of Rayleigh- 
Benard convection rotated about a vertical axis, Cross, et aB computed the local orientation of 
rolls to characterize domain structure. Ouyang and Swinneyil also used the local orientation to 
analyze patterns in a chemical reaction-diffusion system. Gunaratne, et al3 have suggested an 
invariant measure of the disorder of locally striped patterns, and they have demonstrated its use on 
reaction-diffusion patterns. In a set of recent papers, Newell and coworkerst3il3 have begun to use 
local wavenumbers to study the behavior of phase-diffusion equations in the presence of defects. 

Although these previous studies clearly demonstrated the potential utility of local pattern prop- 
erties, they were somewhat limited by the large amount of time it took to process each snapshot 
of a complicated time-dependent behavior. The researchers used a variety of methods including 
nonlinear least squares fits to small patches of the pattern)*!! automatedQ and semi-automatedH de- 
termination of pattern "skeletons" , and modern wavelet-based procedures. EJtl'E^rtlJ In this Letter we 
describe a new, fast method for calculating these local properties based on ratios of simple partial 
derivatives. 

For patterns that are locally striped, we can approximate each field point u{x) using: 

u{x) = A(x) cos(4>(x)). (1) 

The local wavevector k is defined as: 

k(x) = P(p(x). (2) 



Sufficiently far from defects and grain boundaries, we expect that the variations in A(x) are small 
compared to the variations in 4>(x). In that case, the components of the wavevector k can be found 
using simple partial derivatives. For example, 

k l = -^M (3) 

where k x = k • x and d 2 = d 2 /dx 2 . A similar expression can be written for k 2 However, two 
problems are immediately apparent. The first is that when u(x) is close to zero, Eq. (||) will be 
very sensitive to small uncertainties in u{x) due to experimental or numerical noise. This problem 
is easily remedied by taking the ratio of the third partial derivative to the first partial derivative 
for points where u(x) is small. 

The second problem is more subtle. Eq. (^) and the equivalent equation for k 2 yield only the 
magnitudes of the wavevector components k x and k y . From these magnitudes, the vector can be 
specified in any of four directions. We note, however, that periodicity in a real field is invariant 
under rotations by n, so the wavevector (k x , k y ) is equivalent to the wavevector (— k x , —k y ); i.e., the 
wavevector is actually a director, so we only need to determine the relative signs of k x and k y . We 
accomplish this by using combinations of partial derivatives in the x-direction and the y-direction. 
For example, if Eq. (y) is used to determine | A;^ | , the value of k y with its sign relative to k x is 
computed using: 

y d£u(x) 

If d 2 u(x) is close to zero, which occurs when the wavedirector lies in a direction close to the y-axis, 



k y is determined first and then k x . 

Our procedure for determining the values of k x and k y at each point of the field u(x) is sum- 
marized as follows. The global mean of the field, < u(x) >g, is first subtracted from the field. This 
result is then normalized so that the largest absolute value of the field is 1. For each point x, the 
value of the normalized field u(x) and its partial derivatives d x u(x), d y u(x), d 2 u(x), and d 2 u(x) 
are used to choose between the following four cases: 



1. \u(x)\ > max(\d x u(x)\,\d y u(x)\) and \d 2 u(x)\ > \d 2 u(x)\ 

fox = 



d 2 u(x) 



u(x) 
d xy u(x) 



(5) 



ky ~ kx d 2 x u{x) (6) 



2. \u(x)\ > max.(\d x u(x)\,\d y u(x)\) and |<9 2 u(x)| < \dyU(x)\: 
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3. \u(x)\ < mayi(\d x u(x)\,\dyu(x)\) and \d x u(x)\ > \d y u(x)\: 

K x = 



dlu{x) 



k,. 



n'.l; 
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d x u(x) 



0) 
(10) 



4. \u(x)\ < max(\d x u(x)\,\dyu(x)\) and \d x u(x)\ < \d y u(x)\: 
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' dyu(x) 



(11) 

(12) 



If the value of k x is negative, the signs of both k x and k y are changed to ensure that k x > 0, while 
k y contains the information about the orientation of the wavedirector. 

Because this method is local, the effects of noise, of higher harmonics, and of amplitude varia- 
tions neglected in Eq. (|l]) are often noticeable. We reduce these effects by smoothing the k x and 
k y fields over small regions. Although more sophisticated smoothing methods are possible, we have 
successfully employed a simple method in which the k x (or k y ) value is replaced by the average of 
the k x (or k y ) values within a small square region of the field points. We have generally found it 
sufficient to smooth over a region of size A 2 , where A is approximately the average wavelength in 
the pattern. The particular choice of A 2 removes most of the contributions from higher harmonics 
and amplitude variations. Further smoothing can be accomplished by repeating the procedure 
multiple times as necessary. The size of the smoothing region and the number of smoothings can 
be easily determined empirically since local higher harmonics will appear as obvious modulations 
to the wavenumber field. 

As we noted above, the wavevector field is actually a director field, so care must be taken when 
averaging the neighboring k x (k y ) values since k is equivalent to —k. To determine whether to 
include k or —k in the average, we take the dot product of the wavevector at the point of interest, 



ko, and the wavevector of the point to be included in the average, k. If ko ■ k > 0, the angle between 
ko and k is less than n/2 so we average in k x (k y ); however, if ko ■ k < 0, we average in — k x (—k y ). 

Since the local wavenumber is undefined at defects and grain boundaries, our procedure yields 
spurious values of k x and k y at these points. In the averaging procedure we are careful to exclude 
values of k x and k y for which k 2 = k\ + k 2 is far outside the range of expected values. We have 
found that typically only values within a distance of about 1-2 field points around the defects need 
to be excluded from the average. 

Two other factors need to be considered when analyzing experimentally generated data fields (as 
opposed to fields generated by simulating partial differential equations). Typically, experimental 
data have a larger component of large- wavenumber noise than simulational data. The effect of this 
noise can be lessened by applying a Fourier filter to the entire data field before our procedure is 
employed. This Fourier filter should only be applied at wavenumbers well outside the expected 
range of local wavenumbers. Another potential problem is that in some experiments the mean of 
the picture varies spatially, perhaps due to optical imaging effects. One solution to this problem 
is to first determine the mean of the field everywhere through averages over large regions and 
subtract this spatially varying mean from the data field. However, for fields in which the variation 
of the mean is slow compared to typical local wavelengths, a simpler solution is to replace the 
ratios d 2 u(x)/u(x) and dyu(x)/u(x) in Eqs. <^) and (Q) with d^.u(x) / d 2 u(x) and dyu(x) / dyu{x) , 
respectively. 



A final consideration in this procedure is the best method for computing the partial derivatives 
in Eqs. (g)-(12). To obtain high-order derivatives, we perform the differentiation in Fourier space. 



Alternatively, the partial derivatives could be obtained in coordinate space directly using high-order 
finite-difference schemes. We note that the calculation of the derivatives only consumes a small 
fraction of the total time to evaluate a pattern. The computation is heavily dominated by the 
smoothing procedure and by the choice of what derivatives to use at each field point. If, for a given 
application (such as locating disclinations), it is not necessary to determine the wavevector field 
very precisely, little or no smoothing is required and the computation time is reduced drastically. 

Once the values of k x and k y are determined for each point of a data field, local quantities such 
as the wavenumber, the orientation of the periodicity, and the curvature can be computed. As a 
sample of the use of our procedure, we present four figures showing these quantities for patterns 



generated by simulations of the Boussinesq equations describing Rayleigh-Benard convection 



§ 



The local orientation of the rolls, 9{x) = &rctan(k y (x) / k x (x)) , at each point of the pattern 
in Fig. H is shown in Fig. |2[ This figure emphasizes the large regions of straight rolls often found 
within the spiral defect chaos state of Rayleigh-Benard convection. Fig. ^ shows the local wavevector 



magnitude, |fc(x)| = Jk%{x) + ky(x), at each point of Fig. |]. The wavenumber is approximately 
constant across the pattern with a mid-range value (green), but small localized regions of high 
(red) and low (blue) wavenumber are interspersed throughout. (These data, plus the value of the 
roll curvature at each point, were obtained in less than 30 seconds on a single RS/6000 POWER2 
processor running at 66.7 MHz, and even less time would be required at lower levels of smoothing.) 

The region inside the small rectangle in Fig. [T| is shown in greater detail in Fig. |](a). Fig. ||(b) 
and (c) show the local roll orientation and local wavevector magnitude, respectively, for this pattern. 
The rotation of the roll orientation about the centers of each of the two spirals is shown clearly in 
Fig. Q(b). Also, along the left side and along parts of the top and bottom portions of the figure 
are regions of almost constant color indicating that the rolls are approximately straight in those 
regions. In Fig. ||(c), areas with compressed rolls (high wavenumber) are highlighted in red (e.g., 
the bottom left corner), while defects and regions with low wavenumber are shown as blue. 

As an illustration of how these local pattern properties relate to the features of the pattern, 
Fig. ||(a) shows a simple pattern with a single dislocation moving up the page. Fig. ||(b) an d (c) 
show the magnitude of the wavevector and roll curvature C = V • k, respectively, as a function of 
position for this pattern. The deformation of the pattern caused by the dislocation is clear from 
this figure. A "wake" of locally high wavenumber follows the defect while small regions of low 
wavenumber form diagonally in front of the defect. The locality of the majority of the deformation 
is also evidenced by the small regions of high curvature shown in Fig. ||(c). 

We are currently studying the use of statistics based on the quantities described above as order 
parameters for describing transitions between different spatiotemporal chaotic states. In addition, 
we are looking at how the spatial and temporal behavior of the local quantities are related to the 
dynamics of the chaotic states. Of particular interest is the behavior of defects and disclinations, 
which can be located and classified based on the local wavevector field. These works will be 
described in future publications. 

We have presented a novel, fast, automated method for determining local pattern properties as 



a function of spatial position for patterns that are locally striped. Although we have emphasized 
our method's utility for characterizing spatiotemporal chaotic states, we stress that our method is 
simple and general, and it can be applied to locally striped patterns independent of their origin. 
Because of its speed, this new method may allow for the calculation of new statistics to characterize 
disordered patterns in a wide variety of scientific disciplines. 

Correspondence should be directed to egolf@tc.cornell.edu. 
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Figures 

FIG. 1. The temperature field at the midplane of a simulation of the Boussinesq equations with 
reduced Rayleigh number e = 0.6 and Prandtl number a = 1.0 in a box with square periodic 
lateral boundaries with aspect ratio T = 160 x 160 and 512 grid points in each lateral direction. 
Dark points indicate cold downflow; white points indicate warm upflow. 

FIG. 2. The local roll orientation 0(x) = &rctan(k y (x)/k x (x)) in radians for the field shown in 
Fig. |l[ Two smoothings, each with a radius of 3 field points (~ A/2), are performed. 

FIG. 3. The local wavevector magnitude \k(x)\/k c for the field shown in Fig. [j], where k c is the 
magnitude of the wavevector at the onset of convection. Two smoothings, each with a radius 
of 3 field points (~ A/2), are performed. 

FIG. 4. (a)A detailed view of the pattern inside the small rectangle in Fig. [j], (b) its local roll 
orientation, and (c) its local wavevector magnitude \k(x)\/k c . The color schemes for (b) and 
(c) are the same as in Fig. || and Fig. |3|, respectively. 

FIG. 5. (a) A small region (r = 55 x 48) from the temperature field at the midplane of a simulation 
of the Boussinesq equations with reduced Rayleigh number e = 0.55 and Prandtl number a = 1.1 
in a box with square periodic lateral boundaries with aspect ratio T = 201 x 201 and 512 grid 
points in each lateral direction. Dark points indicate cold downflow; white points indicate warm 
upflow. (b) The local wavevector magnitude \k(x)\/k c for the pattern shown in (a). The colors 
range from blue (1.32) to red (1.49). (c) The local roll curvature C = V • k for the pattern 
shown in (a). The colors range from blue (0.00 x 10~ 3 k c ) to red (7.33 x 10 -3 k c ). 
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This figure "figurel.gif" is available in "gif ' format from: 



http://arXiv.org/ps/comp-gas/9702001v2 



This figure "figure2.gif" is available in "gif ' format from: 



http://arXiv.org/ps/comp-gas/9702001v2 



This figure "figure3.gif" is available in "gif ' format from: 



http://arXiv.org/ps/comp-gas/9702001v2 



This figure "figure4.gif" is available in "gif ' format from: 



http://arXiv.org/ps/comp-gas/9702001v2 



This figure "figure5.gif" is available in "gif ' format from: 



http://arXiv.org/ps/comp-gas/9702001v2 



